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Abstract In order to understand the role of space in ecological communities where each 
species produces a certain type of resource and has varying abilities to exploit the resources 
produced by its own species and by the other species, we carry out a comparative study of an 
interacting particle system and its mean-field approximation. For a wide range of parameter 
. values, we show both analytically and numerically that the spatial model results in predictions 

that significantly differ from its nonspatial counterpart, indicating that the use of the mean- 
^ field approach to describe the evolution of communities in which individuals only interact 

D . locally is invalid. In two-species communities, the disagreements between the models appear 

pL^ ' either when both species compete by producing resources that are more beneficial for their 

own species or when both species cooperate by producing resources that are more beneficial for 
the other species. In particular, while both species coexist if and only if they cooperate in the 
mean-field approximation, the inclusion of space in the form of local interactions may prevent 
coexistence even in cooperative communities. Introducing additional species, cooperation is 
no longer the only mechanism that promotes coexistence. We prove that, in three-species 
I communities, coexistence results either from a global cooperative behavior, or from rock-paper- 

scissors type interactions, or from a mixture of these dynamics, which excludes in particular all 
cases in which two species compete. Finally, and more importantly, we show numerically that 
the inclusion of space has antagonistic effects on coexistence depending on the mechanism 
involved, preventing coexistence in the presence of cooperation but promoting coexistence 
in the presence of rock-paper-scissors interactions. Although these results are partly proved 
analytically for both models, we also provide somewhat more explicit heuristic arguments to 
. explain the reason why the models result in different predictions. 

m 

00 

fSJ ' 1. Introduction 

(N ■ 

fSJ I The main objective of tliis article is to understand the role of space, taking the form of local but 

' also stochastic interactions, in the long-term behavior of ecological communities. This is carried 

out through the comparison of analytical and numerical results for an interacting particle system 
and its nonspatial mean-field approximation. Both the spatial and nonspatial models mimic the 
dynamics of communities in which each species produces a certain type of resource and has varying 
r> ■ abilities to exploit the resources produced by its own species and by the other species. As explained 

, later in this paper after the rigorous mathematical description of the interacting particle system 

and its mean-field approximation, this framework can be used to model a wide variety of biological 
interactions such as competition, mutualism, allelopathy, or predation. In the interacting particle 
system, individuals are located on an infinite grid and only have access to the resources produced 
by their nearest neighbors. In contrast, the mean-field approximation assumes that all individuals 
interact globally therefore the dynamics only depend on the overall densities of species, which results 
in a nonspatial deterministic model that consists of a system of coupled differential equations. The 
reason for this comparative study is mainly motivated by the following key question raised in 
the seminal article of Durrett and Levin [5] about the importance of space: which details should 
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be included in a mathematical model and which ones can be ignored? The simplest approach in 
modeling the evolution of inherently spatial interacting populations is to assume that the system is 
homogeneously mixing and to use a system of ordinary differential equations, which is a common 
approach in the life science literature. At the other extreme, one can use the framework of interacting 
particle systems in which individuals are discrete and space is treated explicitly. For a wide range 
of parameter values, we show that both approaches result in different predictions as for the sets 
of species that survive or coexist in the long run. This indicates that, at least for the type of 
interactions considered in this article, the use of the mean-field approach to model the evolution 
of communities in which individuals are static and interact only locally is invalid. In other words, 
local interactions cannot be ignored and should be included in the model. 

Using the terminology of game theory, a species that has a higher ability, respectively a lower 
ability, than all other species to exploit the resources it produces can be seen as a defector, re- 
spectively a cooperator. In the presence of only two species, the spatial and nonspatial models 
strongly disagree when both species are defectors or both species are cooperators. In particular, 
while mutual cooperation always leads to coexistence in the mean-field model, the inclusion of local 
interactions translates into a reduction of the coexistence region. Past research has revealed that 
the inclusion of space can indeed prevent coexistence. This has been proved analytically for a spa- 
tially explicit version of the Lotka-Volterra model introduced in [15] and the non Mendelian diploid 
model introduced in [12]. These two references and the companion paper |llj . which provides rig- 
orous proofs of the analytical results stated in this article for the interacting particle system, focus 
on spatial and nonspatial models with only two types. Our approach is somewhat different. First, 
one of our main objectives is to give more explicit intuitive explanations of the reason why the 
inclusion of local interactions can affect drastically the long-term behavior of ecological communi- 
ties. In particular, heuristic arguments are provided to this extent. Second, the emphasis here is 
on the three-type models. Not only the analysis of the models in the presence of three species is 
somewhat more challenging but also, and more importantly, both the spatial and the nonspatial 
three-type models exhibit new interesting behaviors that cannot be observed when only two species 
are involved. For instance, whereas cooperation is the only mechanism that promotes coexistence in 
the presence of two species, we prove that there are other such mechanisms when additional species 
are introduced. Interestingly, our analysis also indicates that the inclusion of local interactions has 
antagonistic effects on coexistence depending on the mechanism involved, promoting coexistence 
in some situations but preventing coexistence in other situations. In particular, we believe that the 
three-type spatial and nonspatial models can capture most of the important aspects of communities 
involving a large number of species, whereas the two-type models do not. 

2. Models description 

In order to understand the role of local interactions (and stochasticity) in the dynamics of ecological 
communities including varying abilities to exploit resources, we introduce and investigate a stochas- 
tic spatial model and its deterministic nonspatial analog, and then confront the predictions based 
on both models. The first model is an example of interacting particle system while the second one is 
its mean-field approximation which consists of a system of differential equations. For a description 
of the theoretical framework of interacting particle systems along with biological motivations, we 
refer the reader to Durrett and Levin [4J and Neuhauser |14] . For the sake of concreteness, we think 
of both models as describing the dynamics of ecological communities with n species that we label 
from species 1 to species n, but we point out that these models can also be seen as more general 
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models in game theory with potential applications in population genetics, economics, or political 
sciences. In the ecological context, the individuals of species j has to be thought of as producing a 
resource that we simply call resource j. The interactions, either local in the spatial model or global 
in the nonspatial one, are dictated by the n x n matrix 



M 



( ai,2 ••• Ol,n\ 

02,1 Cl2,2 • • • 0-2,n 



where coefficient ajj- represents the ability of an individual of species i to exploit resource j. It is 
assumed that matrix M has only nonnegative entries and, to avoid degenerated cases later, that 
each column of the matrix has at least one positive entry. Wc note that each row can be seen as 
the abilities of a given species to exploit each of the n resources, and each column as the abilities 
of each of the n species to exploit a given resource. In particular, the assumption that each column 
of the matrix has at least one positive entry means from an ecological point of view that each of 
the resources can be exploited by at least one species. Since different species have a priori different 
abilities to exploit each of the resources, their fitness strongly depends upon the configuration of 
their environment as well. Therefore, there is a constant feedback between the set of fitnesses and 
the configuration of the environment, which creates nontrivial dynamics. 

The spatial model. Following the traditional framework of interacting particle systems, we 
assume that space is discrete and time is continuous. More precisely, the spatial structure is repre- 
sented by the d-dimensional integer lattice, Z'^, while the temporal structure is represented by the 
variable t which is any nonnegative real number. Each lattice site has to be thought of as a spatial 
location which, at any time, is occupied by an individual of one of the n species. This corresponds 
to a high-density limit in which, at death, an individual is instantaneously replaced by the offspring 
of a nearby individual. Hence, the spatial configuration of the system at time t is represented by a 
function r\t that maps the integer lattice into the set of species, with r\t{x) indicating the species of 
the individual at site x at time t. 

To describe the dynamics of the spatial model but also to derive its mean-field approximation, it 
is convenient to think of the spatial structure as a graph. In the spatial model, each lattice site is 
connected by an edge to each of its 2d nearest neighbors, therefore, writing x ~ y to indicate that 
sites X and y are nearest neighbors, the integer 

Ni(x) = card {y G Z*^ : y ~ x and rit{y) = i} 

denotes the number of type i neighbors of site x. The species at each site is updated at the arrival 
times of a Poisson process with intensity 1, that is the times between consecutive updates at a 
given site are independent and exponentially distributed random variables with parameter 1. The 
dynamics are dictated by stochastic and local interactions, meaning that the new type at each 
update is chosen randomly from the neighborhood, in order to model the inclusion of an explicit 
space. More precisely, given that the species at site x is of type j at the time of an update, the new 
species at this site is chosen to be species i with probability 

aijNiix) 



aij Ni{x) + a2,j N2{x)-\ h a„j iV„(x) 
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and we assume, in the event that the denominator is equal to zero, that the update is canceled 
so the species at site x remains unchanged. This probability can be interpreted as follows. Since 
the coefficient aij represents the ability of an individual of species i to exploit resource j, the 
numerator is simply the overall ability of the type i neighbors to exploit the resource at site x at 
the time of the update. Similarly, the denominator represents the overall ability of all the neighbors 
to exploit the resource at site x. Therefore, the probability in ([T]) is the relative ability of the type 
i neighbors to exploit the resource at site x, which we also naturally consider as the probability of 
the new species at site x to be species i. Hence, the new species at each update is chosen randomly 
from the local neighborhood with a selective advantage for neighbors that have a higher ability to 
exploit the resource present at the site to be updated. Finally, we point out that, because there 
are countably many lattice sites, updates at different sites cannot occur simultaneously, so the new 
type can indeed be chosen from the local neighborhood without ambiguity. 

The mean-field approximation. The mean-field model is naturally derived from the original 
interacting particle system by excluding both space and stochasticity. For a discussion about the 
connections between interacting particle systems and their mean-field approximations, we refer the 
reader to Durrett and Levin [5j. To motivate the definition of the mean- field model, we first replace 
the infinite lattice by a finite complete graph, that is a finite graph in which each site is connected 
to all other sites including itself, so that all sites have the same neighborhood that consists of all the 
vertex set. Describing the interactions as previously but using this universal neighborhood rather 
than local neighborhoods results in a nonspatial environment in which the types at different sites 
are independent and identically distributed random variables. In particular, 

Ui{t) = P {rit{x) = i) for all t > and i = 1, 2, . . . , n 

where x is a denominated site, are probabilities that indeed no longer depend on x. The mean- 
field approximation is then defined as the deterministic system of coupled ordinary differential 
equations that describes the evolution of these probabilities. To derive this deterministic system 
from the stochastic dynamics, we first observe that 



Ui{t + h) - Ui{t) = P {7]t+h{x) = i, r]t{x) ^ i) - P {r]t+h{x) ^ i, r]t{x) = i) 

= J2 ^ i'^t+hi^) = ^' '?*(^) ^ (^t+hix) = j, Vt{x) 



(2) 



Moreover, the probability that, at the next update, a type j is replaced by a type i is also the 
probability Uj that the site chosen to be updated is of type j times the probability ([T]) that the new 
type selected is type i. Since in addition sites are independent and the probability of an update at 
site X in a very short time interval is roughly equal to the length of this interval, we obtain 



P iVt+hix) = i, mix) =j) = h u,{t) X + oih) (3) 



ai,j Ui{t) 
aijui{t) H h an,j Unit) 

for /i > small. Substituting ([3]) in (j2]), dividing by h, and taking the limit as /i — >■ 0, we get 



^ = y '^hl^ X n,- - y X (4) 

dt ^ ai j ni H h an,j Un ^ ai,i ui H h a„,i Un 



for all i = 1,2, ... ,n. Each fraction in the first sum of ^ is the relative ability of species i to 
exploit resource j over all the community which, when different species have different abilities to 
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exploit each of the resources, strongly depends on the species distribution. Therefore, the first sum 
can be seen as the overall birth rate of individuals of species i. Note also that the second sum is 
simply equal to Ui, so the mean-field approximation reduces to the system 

*^^(± ^ l)n fori=l,2,...,n. (5) 

The minus one term indicates that the per capita death rate of each species is equal to 1, which 
results from the fact that, in the original stochastic process, each lattice site is updated at the 
arrival times of a Poisson process with intensity 1. 

The importance of space. Even though, following the traditional terminology, the mean-field 
model is called nonspatial deterministic approximation, both models, in fact, include space, but 
at different levels of details. In the interacting particle system, the state space consists of spatial 
configurations of species on the lattice that can only interact with their nearest neighbors therefore 
space takes the form of local interactions, which is referred to as explicit space. The mean-field model 
is obtained by replacing the lattice by a complete graph, in particular space is again included but 
it now takes the form of global interactions, which is referred to as implicit space. The interacting 
particle system describes more suitably populations in which individuals are static and can only 
disperse their offspring over short distances, whereas the mean-field model is more appropriate to 
mimic populations in which individuals either move or have the ability to disperse their offspring 
over very long distances. Most mathematical models introduced in the life science literature that 
describe spatially explicit interacting populations consist however of systems of ordinary differential 
equations that assume that populations are well-mixed. The objective of this study is to prove both 
analytically and numerically that, at least for the type of interactions considered in this article, 
the use of the mean-field approach is invalid for a wide range of parameters in the sense that both 
models result in different predictions. Before stating our results, it is necessary to specify what is 
meant here by: both models result in different predictions. In the interacting particle system, since 
the types at different lattice sites are not independent, some quantities of interest are the so-called 
spatial correlations, the correlations between the types at two sites as a function of the distance 
between these two sites. There is obviously no analog of these quantities in the mean-field model 
since it does not include any geometrical structure. Also, the traditional approach to compare both 
models is to compare the density of a given species in the mean-field model with the probability of 
a denominated site being occupied by this species type in the interacting particle system: whenever 
one of these two quantities converges to zero whereas the other one is bounded away from zero 
uniformly in time, we say that both models result in different predictions. In other words, both 
models are said to result in different predictions for a given set of parameters whenever the set of 
species surviving in the long run differ between both models. 

Biological interactions. The limiting behavior of both the interacting particle system and 
its mean-field approximation depends upon two factors: the initial distribution of species, which 
is a spatial configuration on the lattice for the first model and an n-tuple of densities that sums 
up to one for the second model, and more importantly the matrix M that describes the type of 
biological interactions among species. These different types of biological interactions are modeled 
by the relationships among the matrix coefficients, which can be conveniently described by using 
the terminology of game theory. Also, we say that a species is a cooperator, respectively a defector, 
if it has a lower ability, respectively a better ability, than all other species to exploit the resource 
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it produces. That is, in the column corresponding to the resource produced by this species, the 
smahest, respectively the largest, coefficient is the coefficient on the diagonal. While a species can 
be neither a cooperator nor a defector in communities with more than two species, each species is 
either a cooperator or a defector in its relationship with a second species, which induces three types 
of interactions for each pair of species: we call cooperation the interaction between two cooperators, 
competition the interaction between two defectors, and cheating relationship the interaction between 
a cooperator and a defector in which case the defector is called a cheater. To understand the variety 
of biological interactions that can be encoded by the parameter matrix, it is important to point out 
that, even if the individuals of a certain species do not strictly speaking produce resources that are 
beneficial for their species or the other species, their death produces nevertheless a resource available 
to their neighbors: space. In particular, the matrix coefficient Oij can represent more generally the 
ability of species i to exploit space with whatever an individual of species j produces: space with 
resources that are beneficial for the nearby individuals such as nutrients, or at the opposite space 
with biochemicals that have a harmful effect on the growth of other species and reduce their ability 
to exploit the spatial location available in the context of allelopathic interactions. 

Focusing on interactions involving only two species, when none of the species is a producer, the 
interactions are described by the matrix in which both columns are equal. This models competition 
for space with the common coefficient in each row representing the ability of the corresponding 
species to invade space. The species with the better ability to exploit space is a defector while 
the other species is a cooperator. Competition for space is common in plant communities, but also 
occurs among territorial animals such as the grey wolf, Canis lupus, although in this case this is more 
symptomatic of a defensive behavior. In the former context, interactions among plants that spread 
nearby through their rhizome such as Euphorbia are better captured by the interacting particle 
system while the mean-field model describes more suitably the interactions among seed plants 
that can disperse over large distances through anemochory, zoochory or hydrochory such as the 
dandelion. Taraxacum officinale. When both species indeed produce resources and these resources 
are more beneficial for the other species, a symbiotic relationship called mutualism, both species are 
cooperators: the coefficients off diagonal are larger than the coefficients on the diagonal. Mutualism 
is common for instance in terrestrial plants which live in association with mycorrhizal fungi, with the 
plant providing carbon to the fungus and the fungus providing nutrients to the plant. As mentioned 
above, our models are also suited to mimic allelopathic interactions in which a so-called allelopathic 
species produces biochemicals that have detrimental effects on the growth and reproduction of other 
organisms. In this case and when the other species is indeed susceptible to the allelochemicals, the 
allelopathic species is a defector since it has a better ability to invade a spatial location saturated 
with allelochemicals, while the susceptible species can be either a defector or a cooperator depending 
on its relative ability to invade spatial locations void of allelochemicals. Allelopathic interactions 
are common among invasive plants such as the Spotted Knapweed, Centaurea maculosa, that uses 
biochemicals as a defense against other plants and herbivory. In predator-prey systems, the resource 
that the predator exploits is not produced by the prey, it is the prey itself, therefore one needs to 
slightly reinterpret the microscopic rules of the models: a spatial location is not invaded by a nearby 
predator just after the prey dies but just before, with the death of the prey being caused by the 
predator. In this context, the prey is always a cooperator, while the predator is either a defector 
or a cooperator depending on its relative ability to invade empty spatial locations. Note however 
that our models do not reproduce the oscillations typical of certain predator-prey systems such as 
the one designed in the popular Huffaker's mite experiment [9j since it assumes that the predator 
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can survive in the absence of the prey. A similar choice of parameters can also be used to model 
some cases of facultative parasitism in which the parasite can complete its life cycle without being 
associated with the host, provided one identifies the infected host with the parasite itself, which 
is only valid in contexts where infected hosts are sterile. The analysis of the mean-field model in 
this article and the analysis of the interacting particle system in the companion paper [11] give a 
precise picture of the asymptotic properties of systems involving only two species, which includes 
all the biological interactions mentioned in this paragraph. 

For larger ecological communities involving n species, the appropriate relationships among the 
matrix coefficients can be deduced by looking at the type of biological interactions that relates 
any two species. In this case, the analysis of the mean-field model pays a particular attention to 
coexistence and gives sharp results for communities that exhibit rock-paper-scissors dynamics in 
which three species together may coexist whereas, in the absence of any of these three species, the 
other two species cannot: paper covers rock, scissors cut paper and rock smashes scissors. Such 
dynamics have been observed in the common side-blotched lizard, Uta stansburiana, that exhibits 
color polymorphisms [17]. Rather than including three different species, the system consists of male 
individuals of the common side-blotched lizard with three different phenotypes: orange-throated 
males, yellow stripe throated males, and blue-throated males, corresponding to three different 
mating strategies. While evolutionary theory predicts selection of the best fit, all three phenotypes 
coexist in nature which, according to the model of Sinervo and Lively [17], is due to the fact that 
the three phenotypes interact like in the rock-paper-scissors game. Our numerical and anaytical 
results give a valuable insight into this type of dynamics. 



3. Analytical results for the nonspatial deterministic model 

In this section, we collect a number of analytical results for the mean-field model. We start with 
general lemmas that will be applied repeatedly afterwards to understand the limiting behavior 
of the system with two and three species, respectively. The first step is to identify the positive 
invariant sets of the mean-field model. For I C {1, 2, . . . , n}, / 7^ 0, we let 

Si = |(ui,U2, • • • ,Un) G M" : = for i ^ / and J2iei 

Si = ^{ui,U2, . . . ,Un) & Si : Ui > for i £ 

and simply write Si = Sn and Si = Sn when / = {1,2, ... ,n}. From the fact that Uj can be 
interpreted as a probability, it should be intuitively clear that the set Sn is positive invariant for 
the mean-field model. Also, since individuals of either type cannot appear spontaneously in the 
interacting particle system, the same holds in the mean-field model, so each set Si should be 
positive invariant. These statements are made rigorous in the following lemma. 

Lemma 1 For all I C {1, 2, . . . , n}, 1^0, the set Si is positive invariant, i.e., 

(ui(0), . . . , n„(0)) G Si implies that {ui{t), . . . , Un{t)) G Si for all t > 0. 
Proof. This directly follows from the fact that 

= and Ui{t) = whenever Uj(0) = for all i G {1, 2, . . . , n} 
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and by noticing that the second assertion together with the continuity of the trajectories imphes 
that for any species i, we have Uj(t) > at all times t > 0. This completes the proof. □ 

Lemma 2 For all I C {1, 2, . . . , n}, 1^0, the set Si is positive invariant. 

Proof. In view of Lemma [H it suffices to prove that Ui{t) > whenever Uj(0) > 0. Recalling the 
expression of the derivatives, we have 

^ = (y i]n,> -u,. 

V ~^ ^1 ^ ^ '^nj Un J 

Therefore, we have Ui{t) > Uj(0) e^* > whenever Ui{0) > 0. □ 

Our last preliminary result indicates that, if the ability of species i is strictly larger than the 
ability of species j to exploit resources of either type, then species j is driven to extinction when 
starting with a positive density of each species. The intuition behind this result is that the fitness 
of species i, which is a function of the configuration of species, is always strictly larger than the 
fitness of species j regardless of the configuration of the system. 

Lemma 3 Assume that ai^k > ^or all k = 1,2, ... ,n. Then 

lim Uj(t) = whenever u(0) = (ui(0), . . . ,Un(0)) & Sn- 

t— i> CO 

Proof. Since Sn is positive invariant according to Lemma [21 we have 

n 

Um(t) = 1 and Ui{t) G (0, 1) for alH = 1, 2, . . . , n. 

m=l 

Using in addition that each column of M has at least one positive entry, we obtain 



l<m<n \<m<n 
m=l m=l 



Eam,k Um{t) > and min {a^.fe} < am,k Um{t) < max {am,k}- 

m=l 

In particular, defining 
n 

i*,- = > and i j 



(fli.fe — aj^k) Uk , _ . f ai k — aj k 
and L i j = mm < r- 

ai,kUl-\ \-an,kUn ' l<k<n [ maXi<m<„ |Om,fc } 



k=l 

and using that ^ > aj k for all = 1, 2, . . . , n, we get 



dt\ujj Uj ''^ ~ Uj maxi<rn<n {am,k} ' ~ "-■ 



Since Pj > 0, it follows that 



n,(0) 



lim Ui(t) < , , lim exp(— Pj ,t) = 0. 

t^oo ^ - Ui{0) t^oo ^ 



This completes the proof. □ 
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The nonspatial two-type model. Recall that in the presence of n = 2 species, there are 
only three possible types of interactions: cooperation, competition, and cheating. Interestingly, 
the nonspatial model also exhibits exactly three types of long-term behaviors which are directly 
connected to the type of relationship between the species. More precisely, we have the following 
results: cheating implies victory of the cheater, competition implies bistability, and cooperation 
implies coexistence. To prove these results, we first recall that 



dui ( oi 2 a2i 



dt \ai,2Ul + 0-2,2 U2 ai,i Ul + 02,1 U2 

dU2 ( 02,2 ai,2 



dt V Oi 2 ni + 02,2 U2 Oil U\ + 02,1 U2 



Ul U2 

(6) 

Ul U2- 



By invoking the positive invariance of the set 5*2, we obtain U2 = 1 — ui so the model ^ can be 
simply reduced to the one-dimensional ordinary differential equation 

dui f Oi,2 02,l J 

ui(l-ui). (7) 



dt ""1 + ^2,2 (1 - ^ii) ^^1 + ^2,1 (1 — ""1^ 

By setting the right-hand side of ([7]) equal to zero, we find three solutions for ui, namely, the two 
trivial solutions ui = and ui = 1, and the nontrivial solution 

02,1 (ai,2 - 02,2) 

= 7 -( r- (8) 

ai,2 (02,1 - ai,i) + 02,1 (01,2 - 02,2) 

To express the equilibria of ([6]), we also let 

1 - 01,2(02,1-01,1) 
U2 = l-ui = — (9) 

ai,2 (02,1 - ai,ij + 02,1 (01,2 - 02,2) 

and observe that ei,2 = (^1,-02) € (0, 1)^ if and only if 

(«i,i < «2,i and 02,2 < 01,2) or (oi,i > 02,1 and 02,2 > 01,2). (10) 

Note that the first condition in (jlOp indicates that both species cooperate whereas the second 
condition indicates that they compete. In conclusion, the nonspatial two-type model ([6]) always 
has ei = (1,0) and 62 = (0, 1) as its two trivial equilibria. In the presence of either cooperation or 
competition, there is a third nontrivial equilibrium given by ei,2 = (^1,^2). We now investigate the 
global stability of these three equilibria in details. 



Theorem 4 (Global dynamics) For the mean-field model we have 

1. Species 1 is a cheater: if ai,i > 02,1 and 01,2 > 02,2 then ei is globally stable, 62 is unstable, 
and there is no interior equilibrium. 

2. Species 2 is a cheater: if ai,i < 02,1 and 01,2 < 02,2 then ei is unstable, 62 is globally stable, 
and there is no interior equilibrium. 

3. Competition: if ai,i > 02,1 and 02,2 > ai,2 then the interior equilibrium ei,2 is unstable 
whereas the trivial equilibria ei and €2 are locally asymptotically stable, and the system 
converges to e\ whenever ui(0) > ui and to €2 whenever ui(0) < ui. 
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4. Cooperation: if ai,i < 02,1 and 02,2 < «i,2 then the interior equihhrium ei,2 is globally stable 
whereas the trivial equilibria ei and 62 are unstable. 

Proof. The first two statements follow directly from Lemma [3] and the fact that none of the two 
conditions in (jlOp holds. To prove the last two statements, we first observe that (jlOp holds in both 
cases, therefore we have ei^2 £ (0, 1)^ and the global stability of the equilibria can be determined 
by looking at the global stability of their analog for the system ([7]). Now, letting 

h{ui) = ( , r . "^'^ /-I v) ui{l-ui), 

\ai,2Ui + 02,2 (1 - ui) 01,1 ni +a2,i (1 - ui) J 

the stability of ui G {0, are determined by the sign of 



h'(0) = ^^1:111^ h'{l) = ^H:!^^ and 

0-2,1 - ai,i) (ai,2 - 02,2) [oi,2 (ai,i - 02,1) + 02,1 (02,2 - 01,2)] 



h'iui) = 

02,2 - 01,202,1) 

respectively. It is then straightforward to conclude that 

1. When ai.i > 02,1 and 02,2 > ai,2) the trivial equilibria and 1 are locally stable while the 
nontrivial equilibrium ui is unstable. Returning to the mean-field model ([6]) and using that 
the set S2 is positive invariant, we deduce that 

r ( m f^w J^i ^'^^^ ni(0)>'ui 
hm {ui{t),U2{t)) = < /nA / - 

t^oo I 62 when ui[{j) < Ui. 

2. When ai^i < a2,i and 02,2 < ai,2) the trivial equilibria and 1 are unstable while the nontrivial 
equilibrium ui is stable. Returning to the mean-field model ([6]) and using that S2 is positive 
invariant, we deduce that for any initial condition in S2 we have 

lim ui{t) = ui and lim U2{t) = U2- 
This completes the proof of the theorem. □ 



Note that the results of Theorem H] can simply be expressed in terms of the relative ability of 
each species to exploit the resource it produces, defined as 

e, = and 02 = (11) 

01,1 + 02,1 01,2 + 02,2 

since species i is a cooperator whenever 9i < 1/2 and a defector whenever Oi > 1/2. For a summary 
of the results of Theorem HI we refer the reader to the right-hand picture of Figure [3] where the 
phase diagram of the nonspatial two-type model is represented in the ^1-^2 plane. 

The nonspatial three-type model. Note that in the presence of three species, a species can- 
not simply be classified as a cooperator or a defector since it can behave as a cooperator with 
respect to one species but as a defector with respect to another species. This aspect gives rise to 
very rich dynamics. In particular, while cooperation is the only mechanism that leads to coexis- 
tence in communities with two species, there are additional mechanisms that promote coexistence 
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when a third species is introduced. In this subsection, we first look at the local stability of the 
trivial boundary equilibria, and the nontrivial boundary equilibria when they exist. Based on the 
analysis of the local stability of the boundary equilibria, we conclude the subsection with general 
conditions for global coexistence of all three species. These analytical results will be used later as a 
guide to perform numerical simulations of the spatial model in specific parameter regions in order 
to understand the role of space on the global dynamics. For simplicity and to fix the ideas, we 
state some of our results assuming that some species play particular roles in the interactions, but 
additional results can be deduced by replacing (1,2,3) with ((t(1), (t(2), ct(3)) where a is any of 
the six permutations of the group 63. To start our analysis of the nonspatial three- type model, we 
recall that the mean-field model with three species is given by 



dui 
Itt 



FAu) 



ail ui 



0-1,2 U2 



— + 

ai,i ui + 02,1 U2 + 03,1 U3 ai,2 ui + 02,2 U2 + 03,2 ^3 

ai,3 U3 



+ 



ai,3 ui + 02,3 U2 + 03,3 -"3 



1 



(12) 



Ui = Gi{u)ui 



where i = 1,2, 3. Note that, since the set 53 is positive invariant, the mean-field model (jl2p can be 
reduced, by letting ^3 = 1 — ui — U2, to the following two-dimensional system: 



dui 

-dF = ^^("^ 



Oj ,1 Ui 



01,1 ui +a2,iU2 + a3,i (1 - 

ai,2U2 



Ui - U2) 



+ 



ai,2 ^^1 + 02,2 U2 + 03,2 (1 — Ui 
ai,3 (1 - ni - U2) 



U2) 



(13) 



+ 



ai,3 ui + 02,3 U2 + 03,3 (1 - - U2) 



Ui = Gi{u)ui 



where i = 1,2. Note also that the analysis of the two- type nonspatial model in the previous 
subsection directly gives the expression of the boundary equilibria of the three-type nonspatial 
model as well as the conditions under which these equilibria indeed exist. There are always three 
trivial boundary equilibria given by 

ei = (1,0,0) 62 = (0,1,0) and 63 = (0,0,1) 

respectively. By setting the right-hand side of the three equations in (|12p equal to zero, we find 
three additional solutions, namely 



ei,2 = (""1,2, ■"2,1, 0) 62,3 = (0, ■U2,3, ■"3,2) 

where for all i ^ j, the density Ujj- is given by 



and 



63,1 = (^il,3,0,U3,i) 



(oj,- a,-,-) 



According to (fTO|) . ior j = i + 1 mod 3, the solution Cij is indeed a nontrivial boundary equilibrium 
in the sense that both densities Uij and Uj^i lie in the interval (0, 1) if and only if 



(a 



< o,j^i and Qjj < Oj,j) 



or 



(aj,j > Qj^i and ajj > Uij] 
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The next theorem shows the following connection between the three-type and two-type nonspatial 
models: when one species, say species 1, has a strictly higher ability than another species, say 
species 3, to exploit resources of either type, species 3 goes extinct and species 1 and 2 interact as 
described by the two-type model (0), just as in the absence of type 3. 

Theorem 5 Assume that ai^k > fl3,fc for all k = 1,2, 3. Tiien for any initial condition taken in the 
positive invariant set S3 we have the following alternative. 

1. Species 1 is a cheater: if ai^i > 02,1 and ai^2 > 0-2,2 then the trajectory converges to ei. 

2. Species 2 is a cheater: if ai^i < 02,1 and 01^2 < «2,2 then the trajectory converges to 62- 

3. Competition: if ai,i > 02,1 and 02,2 > 01,2 then the trajectory of (fT2]) converges to either ei 
or 62 depending on the initial condition. 

4. Cooperation: if ai^i < a2,i and 02,2 < «i,2 then the trajectory of (fT2]) converges to e\^2- 

Proof. In view of Lemma [3] and the fact that the set S3 is positive invariant, for any e > there 
exists a time T > that only depends on e and the initial condition such that 



The main idea of the proof is to use the fact that after entering the set S3 the three-type system 
behaves similarly to the two-type system ([6j) when the parameter e is sufficiently small, and then 
to apply Theorem HI Since all four cases can be treated similarly, we only provide the details of the 
proof for the first case. By applying the chain rule, we obtain 



Then, we distinguish two cases. 

1. If ai,3 > 02,3 then Lemma [3] implies that X\m.t^ooU2 = 0. 

2. If ai_3 < 02,3 then, using that ai^k > for all /i; = 1, 2, 3, we obtain 

«i,3 ""1(0 + a2,3ii2(t) + 03,3 -0,3 (t) > ai,3iii(t) + 02,3^2 (i) > (1 - e) min{ai,3,a2,3} 
for all times t > T. Therefore, we have 



U3(t) < e and ui{t) + U2{t) > 1 — e for all t > T. 



In other words, any trajectory starting in S3 is trapped after a finite time into the set 



•5*3 = {(ni,ii2,ii3) G [0, 1]^ : Ml + U2 + ii3 = 1 and U3 < e}. 





d / ui 
dt \u2 



) 



> — 




} 



+ 



(1 - e) min {01,3,02,3} 



(ai,3 - 02,3) g 



Sinc6 £ can be arbitrary snicill cuid cii^k — ^2,k for k — 1? 2, W6 havG lim^^QQ U2 = 0. 
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In both cases, the trajectory converges to ei. This completes the proof. □ 



To identify general conditions for tristability on the one hand, and coexistence of all three species 
on the other hand, we now look at the local stability of the trivial boundary equilibria. 

Theorem 6 (Local stability of ei) We have the foUowing alternative: 

1. If oi^i > 02,1 9,nd ai 1 > 03 i then ei is a sink, i.e., ei is locaUy asymptotically stable. 

2. If ai,i < 02,1 and ai,i < 03,1 then ei is a source. 

3. if 02,1 < ai,i < 03,1 or 02,1 > ai,i > 03,1 then ei is a saddle node. 

Proof. The Jacobian matrix of the system ()12p evaluated at equilibrium ei is given by 



/ Q Qi,i - 0-2,1 Qi,i - Q3,i \ 



V 



The eigenvalues of (jl4p are given by 

Ai = A2 

with respective eigenspaces 

^Ai =Span((l,0,0)) 



01,1 
Q2,i - ai,i 
ai,i 





02,1 - ai 1 



ai,i 


Q3,l - Ql,l 



(14) 



ai,i 



Q3,i - 
ail 



Ex, = Span((-1, 1, 0)) Ex, = Span((-1, 0, 1)). 

In particular, for j = 2, 3, species j can invade species 1 in the invariant manifold ui + Uj = 1 
whenever oi^i < aj^i while it cannot whenever ai_i > aj^i. The three statements follow. □ 

From Theorem El one deduces immediately a general condition under which the nonspatial three- 
type system is tristable, which is given by Corollary [7] below. More importantly, this theorem allows 
to identify situations in which there is strong coexistence of all three species in terms of permanence 
of the system, whereas any two species cannot coexist in the absence of the third one. This indicates 
in particular that, in the presence of three species, coexistence is possible even in the absence of 
global cooperation, thus the existence of additional mechanisms that promote coexistence while 
the number of species increases. A more detailed discussion about this important result, which is 
stated below in Theorem [51 is given at the end of this section. 

Corollary 7 (Tristability) All three trivial boundary equilibria ei, 62 and 63 are simultaneously 
locally asymptotically stable whenever Ui^i > aj^i for all i ^ j. 

Proof. This follows directly from Theorem [6l □ 

Numerical simulations further indicate that, under the conditions of Corollary [71 for almost all 
initial conditions, i.e., all initial conditions excluding the ones that belong to a set with Lebesgue 
measure zero, the system (|12p converges to one of the three trivial boundary equilibria. The limiting 
behavior of the system depends on both the parameters and the initial condition. 
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Theorem 8 (Heteroclinic cycle) If the following inequalities hold 

0.3,1 < ai,i < a2,i ai,2 < «2,2 < 03,2 02,3 < 03,3 < ^1,3 

62 63 



(15) 



then there is a heteroclinic cycle connecting e\ 
is locally asymptotically stable whenever 



oi,i > 



02,1 + 03,1 



02,2 > 



01,2 + 030 



ei- Moreover, the heteroclinic cycle 
ai,3 + «2,3 



2 2 
while the heteroclinic cycle is repelling whenever 



03,3 > 



oi,i < 



02,1 + 03,1 



02,2 < 



Ol,2 + 0,30 



03,3 < 



ai,3 + 02,3 



(16) 



(17) 



which implies that the system (jl2p is permanent, i.e., there exists a compact interior attr actor that 
attracts all points in the positive invariant set S3 . 

Proof. By Theorem m the mequahties in (|15p indicate that the system (jl2p only has the three 
trivial boundary equilibria as its boundary equilibria. Moreover, from the proof of Theorem [H we 
know that the transversal eigenvalue at equilibrium Ci in direction uj is given by 



OiA 



In order to obtain a heteroclinic cycle ei ^ 62 — ?• 63 ^ ei, each trivial equilibrium must have one 
positive and one negative transversal eigenvalue, and they must be arranged in a cyclic pattern. In 
particular, (|15p is a sufficient condition for the existence of a simple heteroclinic cycle. To prove 
the second part of the theorem, we let F be a heteroclinic cycle which consists of the three trivial 
boundary equilibria 61,62,63 connected by heteroclinic orbits. Following the approach of Hofbauer 
and Sigmund [7| , we define the characteristic matrix of the heteroclinic cycle F as the matrix whose 
entry in row i and column j is the external eigenvalue Gj{ei). This matrix is given by 



C 





Ol,2 — 02,2 

02,2 
0\,3 — O33 



Q2,i - ox,\ 03,1 - 0\,\ \ 



\ 

Then, defining v = (1, 1, 1), we have 



03,3 



01,1 


"2,3 ~ Q3,3 
03,3 



Ol,l 
03,2 — 02,2 
02,2 





I Gi(6i) + G2(ei) + G3(ei) 

Gl(62)+G2(e2) + 03(62) 
V G'i(63)+G2(63) + 03(63) 



/ 02,1 + 03,1 - 


2ai,i 


\ 


Ol,l 






ai,2 + a3,2 - 


2a2,2 




02,2 






^2,3 + ai,3 - 


2a3,3 




\ 03,3 




/ 



(18) 



By Theorem 17.5.1 in Hofbauer and Sigmund [7], the stability of the heteroclinic cycle can be 
determined based on the sign of the coordinates of the vector in (jlSp . More precisely. 
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1. If all the coordinates of the vector (jlSh are positive, then the heteroclinic cycle T is repelling. 
In this case, the system (jl2|) is permanent. 

2. If all the coordinates of the vector (jlSp are negative, then the heteroclinic cycle T is asymp- 
totically stable in S3. In this case, T is an attr actor of the system (fT2|) . 

This completes the proof. □ 



Under the conditions in (115p . there is a weak form of coexistence of the three species in the sense 
that none of the densities converges to zero. However, the trajectories may get closer and closer 
to the boundaries. Under the additional conditions in ()17p there is a strong form of coexistence in 
the sense that the theorem guarantees that the trajectories are eventually bounded away from the 
boundaries: the limit inferior of each of the three densities is larger than a positive constant. To 
find general conditions under which species coexist as in the two-type model due to a cooperative 
behavior, we now look at the stability of the nontrivial boundary equilibria. 

Theorem 9 (Local stability of ei ^) Assume that 

{ai,i < 0-2,1 and 02,2 < 01,2) or (ai^i > 02,1 and 02,2 > 01,2) (19) 
so that ui^2 G (0, 1) and U2,i € (0, 1) according to Theorein\^ Fix 

Ai,2 = (203,1 - ai,i - a2,i)(ai,2 - 02,2) + (2a3,2 - ai,2 - a2,2)(a2,i - ai,i) 
a quantity that we call the invadihility of species 3. Then, we have the following. 

1. If ai^i < 02,1 and 022 < 01,2 and Ai,2 < 0, then ei,2 is locally asymptotically stable. In 
particular, species 3 cannot invade species 1 and 2 in their equilibrium. 

2. If ai^i < a2,i and 022 < ai,2 and Ai,2 > 0, then ei,2 is repelling. In particular, species 3 can 
invade species 1 and 2 in their equilibrium. 

3. If oi^i > 02,1 and 022 > 01,2 and Ai,2 < 0, then ei,2 is a source. In particular, species 3 can 
invade species 1 and 2 in their equilibrium. 

4. If ai,i > a2,i and 022 > ai,2 and Ai,2 > 0, then ei,2 is a saddle node. In particular, species 3 
cannot invade species 1 and 2 in their equilibrium. 



Proof. Recall that under the first condition in ()19p . the nontrivial boundary equilibrium ei,2 is 
globally stable (cooperation) in the two-dimensional manifold 

{{ui,U2,U3) G M+ : ui,U2 G (0, 1) and ui + U2 = 1} 

whereas under the second condition in (jl9p it is unstable (competition) in this manifold. To have 
a complete picture of the local stability of ei,2 and deduce the invadability of species 3, we look at 
the Jacobian matrix associated with the system ([12]) at point ei,2 



/ OF, 



dui 
dF2 



1 / 

- (ei,2 



(ei,2 



dui 
OF. 
\ 9^ ^''''^ 



dFi 
du2 
dF2 
du2 

du2 



(ei,2) 
(ei,2) 
(ei,2) 



dFi 
du3 
dF2 
du3 
dFs 
du3 



(ei,2) 
(ei,2) 
(ei,2) 
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To express the coefficients of tlie matrix, we introduce 

^3 = iO'i,jUi+a2,jU2 + a3jU3){ei^2) = aij ui^2 + a2,j U2,i forj = 1,2,3. 
Since ei,2 is a boundary equilibrium, Gi(ei,2) = ^2(61,2) = 0. Therefore 



dFi 
dui 

dm 



(ei,2) = 
(ei,2) = 



-^2 



'^2,1 01,2 02,2 



-^2 



fXl 2 ii2,l 



«2,1 



dFi 
du2 

dF2 

duo 



(ei,2) = - 
(ei,2) = - 



Ol,l O2 1 0]_ 2 



-^2 



^1,2 



•^2,1 01,2 02,2 



^2 



Ul,2 U2,l 



Using that the third coordinate of ei_2 is equal to zero, we also have 

5^ ^''''^ = 8^2 ^''''^ = ' 

^^3 . ^ _ r ^ - "3,1 ^1,2 03,2 ^2,1 ^ 

Note that the determinant of the first 2x2 block of the Jacobian matrix J^^ ^ is given by 

5Fi dFo ^F^ dFo 

5;^(e..)x^(e„)-^(e,,)x^(e:,)=0 



while the trace of the matrix is given by 

dFi dF2 
d^, ^''''^ + d^2 ^''''^ = 



02,1(01,1-02,1) 01,2(02,2-01,2) 



+ 



■01,2 ^i2,l 



from which we deduce that the eigenvalues are given by Ai = 
A2 -- 



02,1 (01,1 - 02,1) _^ 01,2 (02,2 - 01,2) 



-^2 



- - J X O31U12 , 03,2^2,1 1 

«i,2 -^2,1 and A3 = ' H - 1 

L>i U2 



The eigenspace associated with the eigenvalue is given by Span(ei,2). Since the second eigenvalue 
only depends on the first 2x2 block of the matrix, we expect that the associated eigenspace is 
spanned by (1,-1,0). To prove this assertion with a minimum of calculation, we observe that 



dFi dF2 dF^ 



dFi , , dF2 , , 



since ^1 + ^2 + ^3 = 0. It follows that 

(ei,2) + ^ (ei,2) = 

"2,1 C''Ol U2,l C''Ol 

In particular, recalling that A2 is the trace of the first 2x2 block. 



01,1 Q2,l 



"1,2 

-^2 



■01,2 + 



^2,1 Oi,2 02,2 



-^2 



U2A 



dFi 
du\ 



i,ei,2) 



dFi 
duo 



(ei,2) 



oi,i 02,1 Oj_ 2 



Ol,l 02,1 _ "'1,2 



U10U2A - 



-^2 



■01,2 



"■2,1 _ Oi,2 02,2 



^^2,1^*1,2 = A2 
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Similarly, we prove that 

dF2 dF2 dFi dF2 

therefore, the eigenspace associated with the eigenvalue A2 is spanned by (1, —1, 0). The sign of this 
eigenvalue thus indicates the stability of ei,2 in the invariant manifold U1+U2 = 1 and as expected 
we find the same conditions as in the case n = 2 given in (jl9p . Rather than computing explicitly 
the third eigenspace, we notice that all three eigenvectors form a basis of the three dimensional 
Euclidean space. This together with the expression of the first two eigenvectors implies that the 
third coordinate of the eigenvector associated with the eigenvalue A3 is different from zero. In 
particular, the sign of A3 indicates whether ei,2 is repelling or not in the direction of 63. To study 
the sign of this eigenvalue, we first observe that, by definition of ei,2, 

02,1 (ai,2^ii,2 + 02,2 ^^2,l) = 01,2(01,1^1,2 + 02,1^^2,1) 
which allows to simplify the expression of A3 and obtain 
sign(A3) = sign [ 2ai,2 03,1 ui,2 + 2a2,i 03,2 ■02,1 - 01,2 (ai,i + 02,1) ^1,2 - 02,1 (ai,2 + 02,2) M2,i ] 
= sign [ai,2 (2a3,i - ai,i - 02,1) ui,2 + 02,1 (203,2 - 01,2 - 02,2) ^^2,1 ] 
= sign [Ai,2 [01,2 (02,1 - ai,i) + 02,1 (01,2 - 02,2)]] • 
The result then follows directly from the previous equation and Theorem HI □ 



Even though the four conditions in Theorem [9] are somewhat complicated to understand due to the 
expression of Ai,2, they can be used to deduce sufficient conditions that are more meaningful from 
an ecological point of view. Since species 3 can invade species 1 and 2 in their equilibrium when 
species 1 and 2 cooperate and Ai,2 > or when species 1 and 2 compete and Ai,2 < 0, we obtain 
that a sufficient condition for invadability is given by 

01,1+02,1 012 + 02,2 
03,1 > ^ and 03,2 > 

In other words, species 3 can invade species 1 and 2 in their equilibrium whenever its ability to 
exploit resource of type 1, respectively 2, is larger than the average ability of the two other species 
to exploit this resource. The last theorem, which is partly based on Theorems [8][9l gives general 
sufficient conditions under which strong coexistence of all three species is possible in the sense that 
the three-type system is permanent. For more convenience, we distinguish four cases depending on 
the number of nontrivial boundary equilibria. We only give particular conditions but, as point out 
above, additional conditions can be deduced by replacing (1, 2, 3) with (cr(l), (t(2), cr(3)) where a is 
any of the six permutations of the group 63. 

Theorem 10 (Coexistence) Let A denote the set of nontrivial boundary equilibria. If conditions 
in one of the following four cases hold, then the system is permanent. 

0. With A = it sufhces that conditions ([E]) and ([IT]) hold. 

1. With A = {61,2} it suffices that Ai,2 > and 



03,1 < ai,i < 02,1 a2,2 < niin(a3,2,ai,2) 02,3 < 03,3 < 01,3 



(20) 
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2. With A = {61,2,62,3} it sufRces that min(Ai,2, A2,3) > and 

^3,1 < < 02,1 02,2 < min(a3,2,ai,2) 03,3 < mm(ai,3, 02,3) (21) 

3. With A = {61,2,62,3,63,1} it suffices that mm(Ai,2, A2,3, A3,i) > and 

ai,i < min(a2,i,a3,i) 02,2 < mm(a3,2, 01,2) 03,3 < min(ai,3, 02,3) (22) 

Proof. Permanence of the system under conditions is simply Theorem [8] that we repeat here 
to have a summary of all the coexistence results. Since permanence can be proved similarly under 
the other three conditions, we only focus on the last case when there are three nontrivial boundary 
equilibria. The proof is based on an application of Theorem 2.5 in Hutson |10] . First, we introduce 
the function P{ui,U2,U3) = ui U2 U3 and observe that this function is equal to zero on the boundary 
of 5*3 and strictly positive in the interior of S3, which makes P an average Lyapunov function of 
the system (jl2|) . From the expression of the derivatives in (|12|) . we then obtain 

, I dP If dui du2 dus \ ^ ^ ^ 

w '■= T, ~nr = ~nr ""2 + -rr u^m + -— uin2 = Gi + G2 + G3. 

P dt Ui U2 U3 V / 

By Theorem |4l the condition (j22p implies that the w-limit sets of the boundary of S3 consists of 
the six trivial and nontrivial boundary equilibria. Therefore, in order to prove that the system 
is permanent, it suffices to prove according to Theorem 2.5 in Hutson [10] that the function 
evaluated at each of the six boundary equilibria is positive. We conclude by observing that 



(Gi + G2 + G3)(6,) = + > 



(Gi+G2 + G3)(6, 



> 



for all i ^ j provided Aj,j > and the conditions in (|22|) hold. □ 

Each of the four sets of conditions in Theorem [TUl which describes the local dynamics at each 
boundary equilibrium, as well as all additional cases that can be deduced from a permutation 
of the three species imply that the system is permanent. However, the theorem does not specify 
whether there are indeed matrices that satisfy these conditions that might be conflicting. To have 
a complete picture of the coexistence region of the three-type nonspatial model, the last step is to 
identify four matrices that satisfy the conditions of Theorem [10] and collect all the configurations 
of local dynamics at the boundary that might also produce permanence of the system. The four 
particular scenarios described in Theorem [10] are illustrated in a schematic manner by the four 
pictures in the first row of Figure [1] The second row gives the solution curves of the system when 
the dynamics are described by each of the following four matrices, and it is straightforward to check 
that these matrices indeed satisfy each of the four sets of conditions of the theorem. 



Mn = 4 1 Ml = 2 M2 = 2 1 M- 
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2 2 2 2 




3 3 3 3 



Figure 1. The four possible scenarios. 

To look for all the possible local stability at the boundary equilibria that might produce permanence 
of the system, we first observe that, whenever two species, say species 1 and 2, cooperate with each 
other but are defectors with respect to species 3, we have 

^3,1 < < 0,2,1 and 03,2 < 02,2 < 0-1,2 

which directly implies that Ai^2 < 0. Therefore, species 3 cannot invade species 1 and 2 in their 
equilibrium. Similarly, when species 1 and 2 compete with each other but are cooperators with 
respect to species 3, we always have the condition Ai^2 < 0. This implies that the first two pictures 
of Figure [2] that describe the local stability of the two trivial boundary equilibria ei and 62, and 
the nontrivial boundary equilibrium ei,2 cannot occur. Using the two previous observations and 
the fact that it is necessary in order to have permanence of the system that none of the boundary 
equilibria is locally stable, simple graphical methods lead to only five possible scenarios, plus again 
the ones that can be deduced from a permutation of the three species. These five scenarios are the 
ones depicted in the first row of Figure [1] and in the last picture of Figure [2] but we claim that the 
boundary dynamics shown in this last picture are conflicting. To establish this result, we observe 
that the local stability at the three trivial boundary equilibria leads to 

02,1 < ai,i < «3,i ai,2 < ^2,2 < 03,2 03,3 < min(ai,3, 02,3). (23) 
To check whether 62,3 can be a saddle node, we compute 

A2,3 = (2ai,2 - 02,2 - 03,2) (a2,3 - 03,3) + (2ai,3 - 02,3 - 03,3) (03,2 - 02,2) 

= (2ai,2 - ^2,2 - 13,2) (^2,3 - 03,3) + (2ai,3 - 2a2,3 + 02,3 - 03,3) (^3,2 - 02,2) 
= 2 (02,3 - a3,3)(ai,2 - ^2,2) + 2 (03,2 - a2,2)(ai,3 - 02,3)- 
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2 2 2 




Figure 2. Impossible pictures. 



In particular, under the conditions ([23]) . 62,3 is a saddle node if and only if A2,3 > which, in view 
of the signs of the first three factors, requires that 01^3 > 02,3. We prove similarly that, for the 
nontrivial boundary equilibrium 63^1 to be a saddle node, it is necessary that 01^3 < 02,3, therefore 
the two equilibria cannot both be saddle nodes. In conclusion, not only the four pictures in the first 
row of Figure [1] are possible but also they are the only ones that lead to coexistence. 

Summary of the results. We conclude this section with a brief overview of the analytical 
results collected for the nonspatial two-type and three-type models. 

The two-type model. In the presence of only two species, a defector always outcompetes a coop- 
erator. More generally, regardless of the size of the community, a species that has a higher ability 
than all other species to exploit resources of either type outcompetes the other species, due to 
the fact that it always has the highest fitness regardless of the global species distribution. When 
both species compete in a defector-defector relationship then the system is bistable indicating that 
coexistence is not possible and that the outcome of the competition strongly depends on the initial 
conditions. In contrast, two cooperators always coexist. Motivated by the properties of the two-type 
model, our analysis of the three-type model focuses on the two extreme cases that also appear later 
to be the most interesting ones: tristability and coexistence. 

Tristability. In communities involving three species, each of the three trivial boundary equilibria is 
locally stable if and only if all three pairs of species are in defector-defector relationship. Numerical 
simulations of the nonspatial model suggest that, as in the case of bistability for the two-type 
model, for almost all initial conditions, one species outcompetes the other two ones. 

Coexistence and cooperation. In contrast, if none of the three pairs of species is in a competition 
relationship, coexistence is possible. In particular, as in the two-type model, coexistence occurs in 
situations in which any two species cooperate. However, we point out that, whereas in the two- type 
model cooperation implies coexistence, global cooperation alone does not imply coexistence of three 
species since additional invadability conditions are required. Even though we omit the details of the 
proof here, it can be easily seen that if the cooperation between two species is significantly stronger 
than the cooperation between any of these two species and the third one, then one of the nontrivial 
boundary equilibrium becomes locally stable: indeed, it suffices to take the matrix with zeros on 
the diagonal, with 01^2 = 02,1 > 2, and with ones elsewhere. Numerical simulations further indicate 
that the two strong cooperators outcompete the third species. Hence, that any two species coexist 
in the absence of the third one does not imply coexistence of all three species. 

Coexistence of cheaters. One of the most interesting aspects of the nonspatial model is that, 
whereas in the presence of two species cooperation is the only mechanism that promotes coexistence. 
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in the presence of three species there are other such mechanisms. In particular, three cheaters can 
coexist. This happens when each species is a defector with respect to a second species but a 
cooperator with respect to the third species, a situation that we shall call rock-paper-scissors type 
relationship to draw another analogy with game theory. When scissors are minoritary, paper beats 
rock and increases its density until the environment becomes suitable enough for scissors to expand 
again, and so on. In particular, that any two species cannot coexist does not imply that global 
coexistence is not possible. In rock-paper-scissors relationships, there is always a weak form of 
coexistence in the sense that none of the species is driven to extinction. However, the densities can 
be alternately arbitrarily low and, as in situations in which species cooperate, strong coexistence 
in the form of permanence of the system requires additional conditions. 

Coexistence and competition. Finally, assuming that two species compete leads to the existence 
of at least one locally stable boundary equilibrium. In other words, whenever two species compete, 
all three species cannot coexist in the sense that the system is not permanent. There are only 
four boundary dynamics, along with the ones deduced from a permutation of the three species, 
that lead to strong coexistence: the four dynamics depicted in the first row of Figure [H Therefore, 
coexistence can only occur when all pairs of species are either in a defector-cooperator relationship 
or cooperator-cooperator relationship, and the coexistence region covers all the cases where the 
number of cooperative pairs ranges from zero to three: regardless of the number of nontrivial 
boundary equilibria, permanence of the system is possible. In conclusion, coexistence occurs either 
in the presence of a global cooperation, or in the presence of a rock-paper-scissors type dynamics, 
or dynamics that consist of a mixture of these two extreme cases. 



4. Analytical results for the spatial stochastic model 



This section collects important analytical results for the two-type spatial model. We focus on the 
meaning of these results and only give an intuitive idea of the mathematical proofs. For a rigorous 
analysis of the two-type spatial model, we refer the reader to Lanchier [11]. As previously, we 
explore the phase diagram in the 9i-92 plane where 6i defined in (jlip measures the relative ability 
of species i to exploit the resource it produces, therefore species i is a cooperator or a defector 
depending on whether the parameter 9i is smaller or larger than one half. 

The one-dimensional two-type process. To understand the one-dimensional model, we first 
initiate the process from the deterministic configuration in which all sites to the left of site 0, 
including site 0, are of type 1, and all the other sites are of type 2. The evolution rules imply 
that, at the times the rightmost site of type 1, respectively, leftmost site of type 2, is updated, the 
probability that the site to be updated switches to the other type is given by 

1-9, and =1-9, 



01,1+02,1 01,2+02,2 

respectively. Since in addition all lattice sites are updated at the same rate, we deduce that the 
interface between both types drifts to the left or to the right depending on whether 9i is smaller 
or larger than 92, respectively. Therefore, when 9i < 92 species 2 wins, which is defined for the 
stochastic process as the fact that any lattice site is eventually of type 2, i.e.. 



hm P {rjtix) =2) = 1 for all x G Z. 



(24) 
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Similarly, when 6i > 62 species 1 wins. Standard techniques further imply that the same conclusion 
holds when starting from a configuration with infinitely many individuals of both types, which 
includes all configurations in which the types at different sites are independent and identically 
distributed and in which both types occur with positive probability. 

To understand the neutral case 9i = 82, we initiate the stochastic process from the configuration 
in which sites are independently of type 1 and of type 2 with the same probability. We note that, 
excluding the case 61 = 02 = 1 in which the process is static, the set of interfaces, defined as 
the random set of points in the dual lattice ^ + ^ that lie between two different types, evolves 
almost like a system of independent annihilating symmetric random walks. More precisely, each 
interface jumps one unit to the left or one unit to the right with equal probability except when two 
interfaces are distance one apart. Moreover, when an interface jumps to another interface, both 
interfaces annihilate. Using that symmetric random walks are recurrent in one dimension, it is 
straightforward to conclude that the set of interfaces goes extinct. This translates into a clustering 
of the process, which is defined as the fact that the probability that two lattice sites have different 
types tends to as time goes to infinity, i.e.. 



The behavior of the one-dimensional process can thus be summarized as follows. 

Theorem 11 Assume that d = 1. Then, 

1. Excluding the static case 9i = 62 = 1, the process clusters. 

2. Assume that 9i ^62- Then type i with 9i = max(0i,^2) wins. 

We refer the reader to [TT], Section 2, for a rigorous proof of these two statements. In contrast 
with the predictions based on the mean-field approximation, Theorem [11] indicates that, except 
in the trivial case when the process is static, coexistence is never possible in one dimension even 
when both species cooperate strongly by providing resources that increase the fitness of the other 
species. Moreover, while the mean-field model exhibits bistability when both species are defectors, 
the spatial model always has a "dominant" type when 61 ^62, that is a type that wins even when 
starting at very low density. This dominant species is always the least cooperative one. The most 
interesting behavior appears in the neutral case when the set of interfaces is roughly a collection 
of symmetric annihilating random walks. Results conjectured by Erdos and Ney [6j and proved 
later by Schwartz [16] show that such a system is site recurrent, meaning that, even though the 
set of interfaces goes extinct due to recurrent annihilations, each site of the dual lattice is visited 
infinitely often by an interface. This indicates that the type at each site alternates infinitely often, 
thus giving the impression of coexistence. However, coexistence is not, strictly speaking, possible 
in the sense that there is no equilibrium in which both types are present: except in the static case, 
there are only two stationary distributions, namely the ones that correspond to the configurations 
in which all sites are of type 1 or all sites are of type 2. 

The two-dimensional two-type process. The analysis of the two-dimensional process is more 
difficult since, due to the geometry of the spatial structure, it cannot be reduced to a simple analysis 
of the interfaces. The idea is to look at key particular cases and then apply standard probabilistic 
techniques such as coupling argument, block construction and perturbation argument to obtain 
additional information about the process. Further insight for the two-dimensional process will be 




for all X, y € Z. 



(25) 
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gained through numerical simulations, however, we point out that our analytical results easily 
extend to any dimension d > 1. Precisely, we focus on the two-dimensional two- type stochastic 
model when the parameters are given by each of the following four matrices 

where the parameter e will be typically chosen small and positive. These four matrices represent 
key points of the phase diagram from which we can deduce, using either monotonicity, symmetry 
or perturbation arguments, the behavior of the process in larger parameter regions. 

First, we assume that the local interactions are described by the matrix M4. This case can be 
seen as a case where species 2 is an extreme defector as the second column indicates that this 
species produces resources that are exploited by its own species only, hence sites of type 2 never 
change their type. In contrast, each individual of type 1 has a positive probability to change its 
type at each update provided it has at least one neighbor of type 2. Thinking of site of type 2 as 
occupied and sites of type 1 as empty, the process becomes a pure birth process so it is obvious that 
type 2 wins. Standard arguments allow to extend the result to a certain topological neighborhood 
of the parameter matrix M4 therefore we obtain the following result. 

Theorem 12 For all < 1 there exists 6 < 1 such that type 2 wins whenever 9 < 62 < 1- 

We refer the reader to Section 4, for more details about the proof. By symmetry, we also 
find a parameter region in which individuals of species 1 outcompete individuals of species 2. 
Note that Theorem 1121 holds even when species 2 starts at a very low density, which shows a new 
disagreement between spatial and nonspatial models: similarly to the one-dimensional process, there 
is a parameter region in which species 2 is the dominant species for the spatial model whereas the 
mean-field model is bistable. 

Assuming now that the dynamics are dictated by the matrix M5 with e = 0, the resources are 
equally shared by both types so the new type at each update is simply chosen uniformly at random 
from the neighborhood. This process is known as the voter model introduced in [Tl[8]. There exists a 
duality relationship between the voter model and coalescing random walks that allows to prove that 
clustering occurs in the sense of (|25p in one and two dimensions whereas there exists a stationary 
distribution in which both types coexist in higher dimensions. Taking e > in the second matrix, 
a similar duality relationship between the process and a certain system of random walks allows 
to show that type 2 wins. This result is expected since in this case species 2 is a cheater: it has 
a higher ability than species 1 to exploit resources of either type. Since this holds for arbitrarily 
small parameter e > and the probability that species 2 outcompetes species 1 is nonincreasing 
with respect to 61 and nondecreasing with respect to 02, we obtain the following result. 

Theorem 13 Assume that 9i < ^ < 92- Then, in any dimension, type 2 wins. 

We refer to Section 3 in [11] for a rigorous proof of this theorem. By symmetry, an analogous result 
can be deduced by exchanging the role of species 1 and 2. In particular. Theorem 1131 shows that, in 
any spatial dimension, if one species is a cooperator and the other species a defector, which results 
in the presence of a cheater, then the cheater is always the dominant species. This agrees with the 
predictions based on the mean-field model. 

The analytical results obtained when the interactions are described by the last two matrices are 
more interesting from a biological perspective. First, we observe that, when the parameters are 
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given by the matrix Mg, individuals of either type change their type at each update whenever there 
is at least one individual of the other type in their neighborhood which, in contrast with the voter 
model, leads to dynamics that somewhat favor the local minority. This corresponds to a case of 
extreme cooperation in which each type only provides resources to the other type. This process is 
known as the threshold voter model for which Liggett [13] has proved that coexistence occurs in 
any dimension d > 1. Relying on perturbation techniques, coexistence can be extended to a certain 
topological neighborhood of the matrix Mq. 

Theorem 14 If d > 1 then there is e > such that coexistence occurs when max(0i, 62) < e. 

See Section 5 in |llj for a proof of this theorem. Note that, while Theorem 1141 shows that increasing 
the spatial dimension allows two species to coexist, it also requires the cooperation to be strong 
enough, whereas results for the mean-field model indicate more generally that cooperation always 
leads to coexistence. 

Interestingly, our next result shows that the inclusion of local interactions indeed translates into 
a reduction of the coexistence region. Before stating this result, we first observe that, when the 
interactions are dictated by the matrix My, the process consists of a mixture of a voter model and 
a threshold voter model: at each update, type 1 is replaced by type 2 whenever it has at least one 
type 2 in its neighborhood whereas type 2 is simply replaced by a type chosen uniformly at random 
from its neighborhood. In this case, there is again a duality relationship between the process and a 
certain system of random walks whose analysis reveals that type 2 wins. Moreover, the techniques 
of the proof allow for the application of a perturbation argument, so the same conclusion holds in 
a certain topological neighborhood of the matrix My which induces the existence of a parameter 
region for which coexistence occurs for the mean-field model but not the spatial model. 

Theorem 15 There is e > such that type 2 wins whenever 61 < e and 62 > ^ ~ ^■ 

The proof of this theorem can be found in [11], Section 6. This result indicates that even though 
two species cooperate, they cannot coexist in space if their relationship has too much asymmetry, 
that is when one species cooperates much more than the other. In such situation, the dominant 
species is again the least cooperative one. 

5. The role of space 

In order to understand the role of local interactions (and stochasticity) , we confront the analytical 
results obtained for the nonspatial and spatial models as well as additional numerical results pro- 
vided in this section for the spatial model. Since numerical simulations of stochastic spatial models 
are somewhat difficult to interpret, we start by devoting few lines to mention important theoreti- 
cal results that justify the approach we have followed to draw general conclusions from particular 
simulation results. Numerical and analytical results are used to first obtain a complete picture of 
the differences between both models when only two species are present. We note however that a 
complete understanding of the stochastic three-type model based on simulations is out of reach due 
to a too large number of parameters. Our approach is to use the results obtained for the nonspatial 
model and the spatial two-type model as a guide to find certain three-dimensional submanifolds of 
the set of parameters along which interesting behaviors emerge. Numerical simulations along these 
manifolds are then performed in order to obtain three-dimensional phase diagrams of the models 
when the set of parameters is restricted to particular subsets. 
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Numerical simulations. In the numerical simulations of tlie spatial model, the population 
evolves on a 400 x 400 lattice with periodic boundary conditions. That is, sites on the bottom row 
are neighbors of those on the top row, while sites on the left edge are neighbors of those on the 
right edge. In addition, the process starts from a uniform product measure, meaning that the types 
at different sites are initially independent and of either type with probability 1/n where n denotes 
the total number of species, so sites are equally likely to be of either type. The process always halts 
after exactly 3.2 x 10^ updates which corresponds to approximately 2000 units of time, which is 
simply obtained by dividing the number of updates by the number of sites and using the fact that 
each site is updated in average once per unit of time. 

The first difficulty in interpreting correctly these spatial simulations is that the stochastic process 
on a finite connected graph always converges to a configuration in which all sites have the same 
type, therefore coexistence is not possible on a 400 x 400 lattice, whereas our analytical results 
show that two types can coexist on the infinite lattice. However, rigorous research in the topic of 
interacting particle systems indicates that, in the presence of a dominant type, this type typically 
invades space linearly in all the directions. In contrast, analytical results notably about the voter 
model on the torus [2] show that, in the absence of a dominant type, the time to fixation of the 
process on large but finite graphs is sufficiently large that numerical simulations indeed reflect a 
transient behavior of the finite system that is symptomatic of the long-term behavior of its infinite 
analog. In particular, it is accurate to define a dominant type as a type which is able to invade our 
finite lattice in less than 2000 units of time while starting from a low density. 

The second difficulty is that the absence of a dominant type as defined above, namely the simul- 
taneous survival of at least two types up to time 2000, is not always symptomatic of coexistence 
for the infinite analog. The reason is that, as proved in [2j for the voter model, in case the infinite 
system clusters, the time to fixation of its finite analog is excessively large, which gives the im- 
pression of coexistence, whereas different types cannot coexist at equilibrium even for the infinite 
system. However, spatial correlations emerge quickly enough so that properties of the configuration 
of the finite system at time 2000 allow to speculate about the infinite system. More precisely, in 
the absence of a dominant type, the dichotomy between clustering and coexistence will be based 
on a quantity that we shall call clustering coefficient and that we define as the percentage of edges 
that connect sites of the same type. Note that the clustering coefficient also gives a good approxi- 
mation of the probability that two nearest neighbors are of the same type. To distinguish between 
clustering and coexistence, we invoke a result proved in p| about the diffusive clustering of the two- 
dimensional voter model, which roughly indicates that the clustering of the process is extremely 
slow and motivates the use of the voter model as a test model. Numerical simulations of the voter 
model, which is obtained from our general model by choosing a matrix in which all the coefficients 
are equal, give the values 86% and 81% for the clustering coefficient of the two-type voter model 
and three-type voter model, respectively, at time 2000. Therefore, in the absence of a dominant 
type, we will assume that the infinite system clusters if the clustering coefficient of its finite analog 
is larger or equal to 86% in the presence of two types and larger or equal to 81% in the presence of 
three types, and that the process coexists otherwise. In conclusion, the long-term behavior of the 
infinite system is related to the clustering coefficient of its finite analog as follows: 

1 . coexistence occurs when the clustering coefficient of the finite system is smaller than that of 
the corresponding voter model, that is either 86% or 81%, 

2. clustering occurs when the clustering coefficient of the finite system is simultaneously larger 
than that of the corresponding voter model and different from 100%, 
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3. existence of a dominant type occurs when the clustering coefficient is equal to 100%, 

and we again point out that past research in the field of interacting particle systems allows to spec- 
ulate with some confidence about the long-term behavior of the infinite system based on numerical 
simulations and the previous classification. 

The two-type models. In the presence of two species, the dynamics only depends on two 
parameters, namely, the relative abilities 6i and 02 of species 1 and 2, respectively, to exploit 
resources of their type. Therefore, the long term-behavior of the spatial and nonspatial models can 
be simply summarized in a two-dimensional phase diagram. The phase diagrams of the spatial and 
nonspatial models are depicted on the left-hand side and right-hand of Figure El The phase diagram 
of the nonspatial model is completely understood analytically whereas the one of the spatial model 
is obtained from a combination of analytical and numerical results. Labels on the parameter regions 
surrounded with dashed lines refer to the theorems of Section [H In both diagrams, the different 
phases are delineated with thick continuous lines, and we distinguish three cases. 

Presence of a cheater. In defector-cooperator relationships, in which case we call the defector a 
cheater, analytical results for both the spatial and the nonspatial models indicate that the cheater 
outcompetes the other species (see Theorems HI and I13p. The intuition behind this result is that, 
since the cheater has a better ability to exploit resources of either type, its fitness is always higher 
than the fitness of the other species regardless of the configuration of species. Therefore, even when 
starting at very low density, the cheater always wins. 

Space and competition. In the presence of competition, which we defined as the relationship 
between two defectors, analytical results for the nonspatial deterministic model indicate that the 
system is bistable (see Theorem H]). That is, the limiting behavior depends on the combination 
of both the parameter values and the initial densities. The parameters being fixed, provided the 
initial density of either species exceeds some critical threshold, this species outcompetes the other 
species. In contrast, analytical results for the spatial model (see Theorem I12p reveal the existence of 
parameter regions in which there is a dominant type that outcompetes the other species regardless 
of the initial configuration provided it starts with a positive density. Numerical results for the spatial 
model further indicate that when 9i ^ 02 there is a dominant type which is uniquely determined 
by the parameters of the system: the dominant type is always the least cooperative species (see 
the upper panel of Tabled]). Additional simulations have been performed to check that this type 
indeed outcompetes the other species even when it initially occupies only 5% of the lattice sites. 
In the neutral case, the system clusters. The spatial correlations become stronger and stronger and 
the boundaries of the clusters sharper and sharper as the common value of the relative ability of 
both species to exploit their own resources increases. Observing in addition that under neutrality 
the number of lattice sites of a given type in the finite system is a Martingale, it can be proved 
that the probability that this type wins is simply equal to its initial density. 

We now give a heuristic argument of the reason why competition translates into bistability in 
the absence of space whereas, excluding the neutral case, the dominant type is uniquely determined 
by the parameters of the system in the presence of a spatial structure. We first observe that if a 
species has a very high ability to exploit the resources it produces but a poor ability to exploit 
the other resources then, in the absence of space, the fitness of this species is roughly proportional 
to its density. Therefore, in the presence of competition, species at low density have a negative 
growth rate and species at high density a positive growth rate, which implies bistability. Including 
space in the form of local interactions drastically modifies the long-term behavior of the system 
due to the fact that individuals can only see their nearest neighbors rather than the whole system. 
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More precisely, since an individual of either type can only appear next to a site already occupied 
by this type (the parent site), spatial correlations build up in such a way that the density of 
individuals of a given type seen by an individual of the same type is in average significantly larger 
in the interacting particle system than in the mean-field approximation. Regardless of the global 
densities, this fraction of homologous neighbors is roughly the same for both species at the interface 
between adjacent clusters. Since in addition the resources produced at a site are only available for 
the nearest neighbors, the global densities become irrelevant in determining the fitness of each 
individual and interfaces expand linearly in favor of the least cooperative species: regardless of its 
initial density, this species is the dominant type. 

Space and cooperation. In the presence of cooperation, analytical results for the nonspatial model 
predict coexistence of both species (see Theorem In contrast, analytical results for the spatial 
model indicate that coexistence is never possible in one dimension (see Theorem II ip . In higher 
dimensions, coexistence becomes possible although there are regions in which coexistence occurs 
for the nonspatial model but not for the spatial model (see Theorems O and [15]) . Numerical 
simulations of the two-dimensional finite system suggest that two cooperative species coexist in the 
neutral case so the coexistence region of the spatial model stretches up to the point that corresponds 
to voter model dynamics. As previously mentioned, analytical results for the voter model show that 
this critical point is included in the coexistence region if and only if the spatial dimension is strictly 
larger than two. There is however a wide range of parameters for which coexistence is not possible 
though both species cooperate (see the lower panel of Tabled]). Therefore, the inclusion of space 




Figure 3. Phase diagrams of the spatial and nonspatial two-type models. 



translates into a reduction of the coexistence region in any spatial dimension. This effect is more 
pronounced in low dimensions and we strongly believe that the coexistence region of the spatial 
model approaches the one of the nonspatial model as the dimension tends to infinity. 

The intuition behind the fact that the inclusion of a spatial structure reduces the ability of species 
to coexist follows the lines of the heuristic argument introduced above. In cooperative mean-field 
systems, a species at very low density turns out to be in a very favorable environment due to the 
high density of the other species and so a large amount of resources to exploit: the smaller the 
density of a species, the larger the fitness of this species, which is a basic mechanism that leads to 
coexistence of both types. In contrast, the inclusion of local interactions translates as previously 
into the presence of spatial correlations that make the density of individuals of a given type seen 
by an individual of the other type significantly smaller in the interacting particle system than in 
the mean-field approximation. Therefore, in the presence of a spatial structure, cooperative species 
cannot fully benefit from the resources produced by the other species simply because these resources 
are geographically out of reach. Under neutrality, both species suffer local interactions equally so 
coexistence is possible, but when the parameters get closer to the region where one species is a 
cheater, individuals of the other species are unable to expand whenever they have at least one 
neighbor of their own type: large clusters of this species tend to shrink quickly while small clusters 
of this species dissolve eventually due to stochasticity. Note also that the number of neighbors 
is smaller, and so the fraction of neighbors of the same type larger, in low dimensions hence the 
harmful effect of space on coexistence is more pronounced in low dimensions. 

The three-type models. The analysis of the mean-field approximation reveals that, in the 
presence of three species, the nonspatial model exhibits a wide variety of regimes. However, an- 
alytical and numerical results for the spatial two-type model indicate that the most interesting 
behaviors in terms of the inclusion of space emerge in parameter regions for which there is either 
bistability or coexistence for the mean-field approximation. Therefore, we mainly focus on the anal- 
ogous cases for the three-type models by investigating the effect of space in regions for which there 
is either tristability or coexistence of all three species when space is absent. We note that, whereas 
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in the presence of two species cooperation is a necessary and sufficient condition for coexistence in 
a nonspatial universe, in the presence of three species, there are additional strategies that promote 
global coexistence. For instance, coexistence occurs for what we shall call rock-paper-scissors dy- 
namics in which any two species cannot coexist but all three species can. Interestingly, numerical 
simulations of the spatial three-type model reveal that the inclusion of space does not have the 
same effect on coexistence depending on the strategies species follow, suggesting that some of the 
strategies that promote coexistence are much more stable than others. 

In the presence of two species, each species is either a defector or a cooperator, which results in 
only three types of relationship: cheating, competition, and cooperation. Introducing an additional 
species increases significantly the number of types of relationships since a species can be a defector 
with respect to a second species but a cooperator with respect to the third one as in rock-paper- 
scissors dynamics. In order to extend to the three-type models results obtained for the two-type 
models, we first assume that any two species have the same ability to exploit resources produced 
by the third one so that each species is either a "global" defector or a "global" cooperator. Under 
this assumption, the dynamics only depend on three parameters being described by matrices that 
belong to the three-dimensional manifold that consists of the matrices of the form 

/ 261 1-02 l-^sX 
Ms =[ 1-01 202 1 - 03 . 

1-02 203 J 

Note that species z is a defector if 0i > 1/3 and a cooperator if 0i < 1/3. Following the same 
approach as for the two-type models and motivated by the analytical results for the nonspatial 
three-type model, we distinguish four cases depending on the number of defectors. 

Presence of a cheater. In case one species is a global defector and the other two species are 
global cooperators, so that the defector is a cheater, analytical results for the nonspatial models 
indicate that the cheater outcompetes the other two species (see Theorem [5]) . Analytical results 
for the spatial two-type model easily extend to prove that the same holds including space in the 
form of local interactions. The intuition is the exact same as for the two-type model. The ability 
of the cheater to exploit resources of either type is larger or equal to the ability of the other two 
species to exploit the same resource. Although the cheater has the same ability of any of the other 
species to exploit the resources produced by the third one, since it has a better ability to exploit the 
resources it produces, its fitness is always strictly larger than that of the other species regardless 
of the configuration of the system. 

Presence of two cheaters. This case occurs when there are two global defectors and one global 
cooperator. That space is present or not, the altruist species is driven to extinction. This can 
be understood based again on the same heuristic argument, which indicates that, regardless of 
the configuration of the system, the cooperator always has a strictly smaller fitness than the two 
defectors. The long-term behavior is therefore similar to that of the two-type models in the presence 
of competition. In the mean-field model, one of the two defectors outcompetes the other one, and 
the limiting densities depend on the combination of the parameter values and the initial densities: 
the system is bistable (see Theorem [5]). In the interacting particle system, the dominant type is 
uniquely determined by the parameter values and is always the least cooperative species. In case 
of a symmetric behavior between both defectors, the community clusters. 

Space and competition. Similarly to the two-type model, the most interesting behaviors emerge 
when all three species compete or all three species cooperate. In the presence of competition, that 
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is when all three species are global defectors, the mean-field model is tristable: the trajectories 
converge to one of the three trivial equilibria in which only one species is present, and the limiting 
densities again depend on both the parameter values and the initial densities. Numerical simulations 
of the spatial model (see the upper panel of Table [5]), indicate that the dominant type is again the 
least cooperative species. In case of equality, additional behaviors can appear: clustering of the 
two species that have the (same) highest relative ability to exploit the resources they produce, or 
clustering of all three species in the neutral case. The disagreement between both models in the 
presence of competition follows from the same argument as for the two-type models. The density 
of individuals of a given type seen by an individual of the same type is in average significantly 
larger in the interacting particle system than in the mean-field approximation in such a way that 
global densities become irrelevant and interfaces expand linearly in favor of the species that has 
the largest relative ability to exploit its own resources, i.e., the least cooperative one. 

Coexistence and cooperation. Recall that in the presence of three species, global cooperation 
alone doe not imply coexistence for the nonspatial model, since additional conditions about the 
nontrivial boundary equilibria are required (see statement 3 in Theorem llOp . However, when the 
dynamics are described by the matrix Mg one has 

Ai,2 = (2a3,i - ai^i - a2,i)(ai,2 - 02,2) + (203,2 - ai,2 - a2,2)(a2,i - ai,i) 
= (l-3ei)(l-302) + (l-3e2)(l-30i) = 2(1 -30i)(l -3^2) 

and similar expressions for A2,3 and A3,i. Therefore, in this special case, coexistence occurs in the 
nonspatial mean- field model if and only if species cooperate. In contrast, in the presence of global 
cooperation, the interacting particle system has no less than seven possible regimes that are all 
stable under small perturbations of the three parameters. The existence of all these regimes follows 
heuristically from the same argument that explains the reduction of the coexistence region for the 
spatial two-type model: the density of individuals of a given type seen by an individual of the other 
type is in average significantly smaller in the interacting particle system than in the mean-field 
approximation, which makes species significantly more cooperative than others unable to survive 
in the presence of space. When cooperation is strong or fair among all three species, then global 
coexistence occurs even in the presence of space. However, when two species are significantly less 
cooperative than the third species, these two species can coexist but outcompete the third one (see 
the middle panel of Table [2]). Similarly, when a single species is significantly less cooperative than 
the other two species, this species is the dominant type: it outcompetes the other two species even 
when starting from a low density (see the bottom panel of Table [2]). 

In conclusion, the phase diagram of the nonspatial model when the set of parameters is restricted 
to matrices of the form Mg can be deduced from Theorems [51 and \T0\ and Corollary [3 The unit cube 
representing the possible values of 9i, 62 and ^3 is divided into eight regions which are delimited by 
the three two-dimensional hyperplanes with equation 0j = 1/3. In the presence of three cooperators, 
coexistence occurs. In the presence of two cooperators and one defector, the defector wins. In the 
presence of one cooperator and two defectors, the system is bistable. Finally, in the presence of 
three defectors, the system is tristable. We refer to Figure[5]for the solution curves of the nonspatial 
model in these four cases. Finally, the three-dimensional phase diagram of the spatial model is 
depicted in Figure HI Excluding the cube with length edge 1/3 at the bottom of the left picture, 
the phase diagram is divided into three regions delimited by three hyperplanes that intersect along 
the segment line, denoted by Ai, that connects (1/3,1/3,1/3) and (1,1,1). Along this segment 
line, all three species cluster, along the hyperplanes two of the species outcompete the third one 
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and cluster, and in the three regions dehmited by the three hyperplanes there is a dominant type 
which is uniquely determined by the parameter values. The right picture shows an enlargement of 
the coexistence region which is included in the cube with length edge 1/3 that corresponds to the 
region where all three pairs of species cooperate. This picture is drawn from the simulation results 
of Table [2j Note that the spatial model exhibits seven regimes in this single cubic region. The first 
row of Figure E] shows three snapshots of the spatial model when parameters are taken along the 
straight line A2. As the common value of d\ and Q2 increases, the system crosses three difi^erent 
regimes: first all three species coexist, then species 1 and 2 coexist but outcompete species 3, and 
finally not only species 3 is driven to extinction but also species 1 and 2 cluster. The second row of 
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Figure 5. Solution curves when the dynamics are described by Mg,- 



the same figure shows three snapshots of the spatial model when parameters are taken along the 
straight line Ai. In this neutral case, cooperation allows all three species to coexist while switching 
to competitive relationships, the system clusters, as shown by the last two snapshots. 

Coexistence of cheaters. Recall that, according to Theorem 1101 and the discussion following the 
theorem, coexistence is possible in the absence of space regardless of the number of pairs of species 
that are in cooperative relationship or in defector-cooperator relationship, whereas coexistence 
cannot occur whenever one pair is competing. After looking at the first extreme case where all 
pairs are in a cooperative relationship, in which case the inclusion of space reduces the coexistence 
region, we now look at the other extreme case where all three pairs are in a defector-cooperator 
relationship that results in rock-paper-scissors dynamics. To understand the role of space in this 
case, we have performed numerical simulations of the spatial three-type model when the local 



The role of space in the exploitation of resources 



33 





01 = 00 


_ Q 


— Q 20 


= 30 


= 40 


= 50 


= 60 


^1 = 70 


0^ = 80 


^1 = 90 




= 0.00 




'iU. ± / .4: J. \ \) 1 ) 


'±0. iO.OO \ V 1 ) 


T: 1 .J. 1 .01 yuo ) 


OU.-LO.OU Ij'Jyj 


8-97 (7^ ^ 


uo.-Ly.zo \i ■'-•) 


•J 1 . J.0.Z0 l_ ' "J J 


'iQ-IQ-91 i7?\\ 
tjy.-Ly.zi 1 u J 


■^8-1Q'91 i7?\\ 
uo.-Ly.z J. 1 u J 




— 0.10 


-1 Q./IQ /'fi71 


■ifi-l Q-dzl /'fiY'l 


'^Q-1 Q dn I'fis'i 

■jy.j.y.'iu l^iJOj 


"rtO.ZU.OO l^'Jyj 


'iJ.ZU.OO {'J^j 


ziQ-9n-9Q ^7^ \ 




u^.zi.zu wo J 


uz.zi.zu wo J 


'^■^-99'99 i7W 


e. 


= 0.20 


oi.lo.4y i^oyj 


OO.ZU.40 i^oyj 


oD.zu.4o i,oyj 


oy.zu.oy [I'jj 


An.nn-or^ l'7r\\ 


40.ZZ.OZ [1 i-j 


4D.Zz.0U I, 'ZJ 


/IQ-99.97 ^7^^ 


4y.Z4.ZD '4J 


c;9-9'J'9^ /7i^^ 
iiL.Lo.Lo \ io] 


0, 


= 0.30 


28:18:52 (70) 


29:21:48 (70) 


34:21:44 (70) 


38:22:39 (71) 


38:23:37 (71) 


40:23:36 (72) 


42:25:32 (73) 


43:27:29 (73) 


44:26:28 (74) 


48:25:25 (75) 


0, 


= 0.40 


25:19:54 (71) 


27:22:49 (71) 


31:23:44 (71) 


33:24:42 (72) 


36:24:38 (72) 


38:25:36 (72) 


39:27:33 (73) 


42:26:31 (74) 


41:28:29 (74) 


44:27:28 (75) 


62 


= 0.50 


24:20:55 (72) 


26:21:52 (72) 


29:24:46 (71) 


31:25:43 (72) 


33:25:41 (73) 


34:26:38 (73) 


36:28:35 (74) 


38:28:32 (74) 


42:28:29 (75) 


40:31:27 (75) 


02 


= 0.60 


22:20:57 (73) 


25:22:51 (73) 


26:24:49 (72) 


28:26:45 (72) 


31:27:40 (73) 


33:28:38 (74) 


34:29:36 (74) 


36:30:33 (74) 


36:30:32 (75) 


39:30:29 (76) 


02 


= 0.70 


20:21:57 (74) 


24:23:52 (73) 


24:27:47 (73) 


28:26:45 (74) 


29:28:41 (73) 


30:29:39 (74) 


32:28:38 (74) 


34:31:34 (75) 


33:32:33 (75) 


35:33:30 (76) 


02 


= 0.80 


19:22:58 (75) 


22:25:52 (73) 


22:27:49 (74) 


25:27:46 (74) 


26:28:44 (74) 


29:31:39 (75) 


30:31:37 (75) 


31:32:35 (76) 


32:33:33 (76) 


35:33:31 (76) 


02 


= 0.90 


18:22:59 (75) 


18:25:56 (75) 


22:26:51 (75) 


23:27:49 (75) 


26:30:43 (74) 


26:31:42 (75) 


28:33:38 (76) 


30:32:36 (76) 


31:34:34 (76) 


31:35:32 (77) 



Densities and clustering coefficient when the interactions are described liy the matrix 7\/c) with ^3 = 0.80 

Table 3 



interactions are described by matrices of the form 

lOi 1\ 

Mg = 1 02 . 

\0 1 03/ 

Recall that, according to Theorem [HI when the dynamics are described by the matrix Mg one 
obtains a heteroclinic cycle by taking 0, < 1 for alH = 1,2,3. Recall also that, under the additional 
conditions (jl6p the heteroclinic cycle is locally asymptotically stable whereas it is repelling, inducing 
permanence of the system, under the conditions (|17p . Numerical simulations of the nonspatial model 
further indicate that there is a unique interior equilibrium. When starting from any condition in Sj, 
excluding this interior equilibrium and under the assumptions ()16p . the trajectories get closer and 
closer to the boundaries of the system indicating that for large times at least one of the densities 
is arbitrarily close to zero, whereas under the assumptions (|17p . the trajectories converge to the 
interior equilibrium. In contrast, numerical simulations of the spatial three-type model indicate 
that taking 0j < 1 for all i = 1,2,3, suffices to obtain coexistence of all three species in space. 
These simulation results are reported in Table [3] where ^3 = 0.80. Note that, even when both 
parameters 9i and 82 are strictly larger than one half, in which case (jl6p is satisfied and so the 
nonspatial model is not permanent, all three densities are bounded away from zero and the cluster 
coefficient is strictly smaller than that of the three-type voter model. In conclusion, whereas the 
inclusion of local interactions reduces the coexistence region when all three pairs of species are in 
a cooperator-cooperator relationship, it increases the coexistence region when on the contrary all 
three pairs of species are in a defector-cooperator relationship. The three snapshots of the spatial 
model in Figure [7] indicate that this also holds when starting with two species at low density. 

We finally explain the reason why local interactions promote coexistence in the case of rock- 
paper-scissors dynamics. The dynamical picture of the spatial model at equilibrium shows that small 
clusters keep forming very quickly but then shrink again immediately: residuals of one species, say 
paper, invade nearby clusters of rocks very quickly resulting in clusters of papers which, when they 
reach a critical size, are invaded and destroyed by nearby residuals of scissors, so the system exhibits 
local traveling waves. Focusing on a reasonably small piece of the lattice, one sees a majority of 
rocks, which are suddenly invaded by a wave of papers that become majoritary, which themselves 
are suddenly invaded by a wave of scissors, and so on. In particular, the dynamics predicted by 




Figure 6. Snapshots of the three-type model at time 200 when starting from a product measure in which sites are 
equally likely to be of cither type and when the local interactions arc described by the matrix Mg. In the first row, 
the parameters are di = 62 and 83 = 0, while in the second row, di = 62 = Os- 




(a) time 40 




(b) time 80 




(c) time 120 



Figure 7. Snapshots of the spatial model starting with 90% of rocks (black) and 5% of papers and scissors (grey 
and white) when the local interactions are described by the matrix Mg with 0i = O2 = O3 = 0.80. 
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the nonspatial model under conditions (|15p -(|16p indeed reflects the dynamics of the spatial model 
viewed in a small region of the lattice. However, due to the presence of local interactions, the global 
densities become irrelevant and small pieces of the lattice that are far from each other are almost 
independent in such a way that the nonspatial dynamics are only symptomatic of what happens in 
small spatial regions. In addition, due to stochasticity, there is a strong asynchrony among nearby 
regions so, even if a region is void in one type, this type can quickly reappear. 

6. Discussion 

The long-term behavior of ecological communities involving competitive defector-defector interac- 
tions, cooperative interactions, and defector-cooperator interactions are difficult to predict mainly 
because the fitness of each individual depends on the environment, which induces a constant feed- 
back between the set fitnesses and the configuration of the community. In situations where one 
species has a higher ability than all the other species to exploit resources of either type, the fitness 
of the first species is always strictly larger than that of the other ones regardless of the configuration 
of the community. Excluding such situations, the interacting particle system and its mean-field ap- 
proximation generally disagree. This is due to the fact that the environment seen by each individual 
strongly differs in the presence and in the absence of a spatial structure. 

Regardless of the size of the community, when all pairs of species are in a competitive relationship, 
the mean-field model predicts multistability of the system: one type eventually takes over and this 
dominant type is determined by the combination of both the parameters of the system and the 
initial densities. In contrast, in the presence of local interactions, the dominant type is uniquely 
determined by the parameters. The reason is that, regardless of the global densities, the fractions 
of species seen from the spatial interface between two types are roughly the same for both types 
which results in traveling waves that expand in favor of the dominant type. 

In the presence of only two species, cooperation is the only mechanism that promotes coexistence: 
in the absence of space, coexistence occurs if and only if cooperation occurs. Introducing additional 
species has two important consequences. First, even in the mean-field approximation, global coop- 
eration alone does not result in the coexistence of all species since a too strong cooperative behavior 
between two species can drive the other species to extinction. In particular, that any two species 
can coexist does not imply that species all together can coexist. Second, there are additional mech- 
anisms that promote coexistence. Rock-paper-scissors dynamics in which all pairs of species are in 
a defector-cooperator relationship is an example of such a mechanism. In particular, that any two 
species cannot coexist in the absence of the other species does not imply that species all together 
cannot coexist. Another important aspect is that, at least in the presence of at most three species, 
whenever two species are in a competitive defector-defector relationship, global coexistence is never 
possible. Therefore, the collection of mechanisms that promote coexistence when three species are 
present consist of global cooperative behavior, rock-paper-scissors type dynamics, and mixture of 
these two extreme cases. 

The effect of local interactions in the presence of either global cooperation or rock-paper-scissors 
dynamics is probably the most interesting aspect of our study: in the presence of mutual coopera- 
tive relationships, the inclusion of local interactions translates into a reduction of the coexistence 
region, whereas in the presence of rock-paper-scissors dynamics, the inclusion of local interactions 
translates into an expansion of the coexistence region. The intuition behind the first statement is 
that the fraction of individuals of a given species seen by an individual of another species is sig- 
nificantly smaller in the interacting particle system than in the mean-field model indicating that. 
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in the presence of local interactions, cooperative species cannot fully benefit from the resources 
produced by the other species. This effect is more pronounced in low spatial dimensions. Finally, 
the intuition behind the second statement is that the behavior predicted by the mean-field model 
in the presence of rock-paper-scissors dynamics is only symptomatic of the local behavior of the 
interacting particle system seen in a fairly small region. Since in addition stocliasticity induces a 
strong asynchrony among disjoint spatial regions, when a species is driven to extinction locally, it 
is quickly reintroduced from nearby regions. 

Acknowledgment. The authors would like to thank an anonymous referee for important comments 
that helped to improve the article. 
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